Analyse time-dependent proteome¶
2. Test burn-in period and steady state¶
We want to simulate 5 min intervals in quasi-steady state, so we cannot start with an empty cell; we have to allow for a burn-in period to reach steady state and discard timecourses of this period.
In [14]:
fig, ax = plt.subplots(size, size, figsize=(16, 16))
plt.suptitle('Selected protein timecourses', fontsize=18, fontweight='bold')
for i in range(size):
for j in range(size):
yseries = results[time]['timecourses'][protnames_plot[i][j]]
xseries = np.arange(0., float(len(yseries)), 1.)
# oddly, sometimes the simulation sometimes 'loses' one data point
length = min(len(xseries), len(yseries))
ax[i, j].plot(xseries[:length], yseries[:length])
ax[i, j].set_xticks(np.arange(min(xseries), max(xseries)+1, 200.0))
ax[i, j].set_title(protnames_plot[i][j], fontdict={'fontsize': 12, 'fontweight': 'bold'})
ax[i, j].tick_params(axis='both', which='major', labelsize=10)
Plot increases and check for saturation of "first derivative":
In [15]:
burnin = len(yseries) - 300 # 300 s simulation, rest is burn-in
fig, ax = plt.subplots(size, size, figsize=(16, 16))
plt.suptitle('First derivatives of selected protein timecourses', fontsize=18, fontweight='bold')
for i in range(size):
for j in range(size):
yseries = results[time]['timecourses'][protnames_plot[i][j]]
dy = np.gradient(yseries)
xseries = np.arange(0., float(len(dy)), 1.)
# oddly, sometimes the simulation sometimes 'loses' one data point
length = min(len(xseries), len(yseries))
ax[i, j].hist(xseries, bins=10, weights=dy)
ax[i, j].axvline(burnin, color='red', linestyle='--')
ax[i, j].set_xticks(np.arange(min(xseries), max(xseries)+1, 200.0))
ax[i, j].set_title(protnames_plot[i][j], fontdict={'fontsize': 12, 'fontweight': 'bold'})
ax[i, j].tick_params(axis='both', which='major', labelsize=10)
3. Plot proteomes at different times¶
Prevent output from scrolling:
In [80]:
%%javascript
IPython.OutputArea.auto_scroll_threshold = 999;
In [81]:
times = np.arange(0, 3600, 600)[:4]
#rows = len(times/columns)
for i, tx in enumerate(times):
#fig, ax = plt.subplots(rows, columns, figsize=(16, 16))
for ty in times[i+1:]:
fig = plt.figure()
plt.plot([1e-1, 1e6-100000], [1e-1, 1e6-100000], color='red', linestyle='-', linewidth=0.5)
scatter_dicts(proteome[tx], proteome[ty], xaxislabel="proteome ({} s)".format(tx+300),
yaxislabel="proteome ({} s)".format(ty+300), fontsize=None)
plt.show()
4. Plot translational impact¶
First we check whether the relevant genes are in our data set.
In [12]:
len(results[time]['timecourses'].keys())
Out[12]:
In [13]:
histones = ['YDR224C', 'YBL002W', 'YDR225W', 'YBL003C', 'YOL012C', 'YBR009C', 'YNL030W', 'YBR010W', 'YNL031C', 'YPL127C']
In [14]:
for histone in histones:
print "Histone {} found in data set: {}".format(histone, histone in results[time]['timecourses'])
In [15]:
TFs = ["YNL027W", "YLR113W", "YMR174W", "YER062C", "YHR084W", "YDL056W", "YLR182W", "YER111C", "YLR118W", "YMR043W"]
In [16]:
for TF in TFs:
print "TF {} found in data set: {}".format(TF, TF in results[time]['timecourses'])
Get mRNAs:
In [17]:
transcriptomes_dict = collections.OrderedDict(sorted(
pkl.load((open('../parameters/transcriptome_time_dependent.p'))).items()))
In [18]:
transcriptomes_dict.keys()
Out[18]:
In [19]:
absolute_transcript_count = {time: sum(transcriptomes_dict[time].values()) for time in sorted(transcriptomes_dict.keys())}
In [20]:
plt.scatter(absolute_transcript_count.keys(), absolute_transcript_count.values(), marker='o')
Out[20]:
TODO: this can be made look nicer for sure.
Generate fold-changes for transcriptome:
In [21]:
time = 30 * 60
In [22]:
print len(transcriptomes_dict[time/60])
print len(transcriptomes_dict[0])
print len(transcriptomes_dict[3000/60])
In [23]:
print transcriptomes_dict.keys()
print len(results[2400]['timecourses'].keys()) # some proteins were not translated (probably no transcripts)
First we try only one pair of times:
In [24]:
time_0 = 0
time_t = 30 * 60 # time always in seconds
In [25]:
def transcriptome_fold_change(time_0, time_t):
tmp = {}
for gene in transcriptomes_dict[time_0/60]:
if gene in transcriptomes_dict[time_t/60] and transcriptomes_dict[time_0/60][gene] != 0 and \
transcriptomes_dict[time_t/60][gene] != 0:
tmp[gene] = np.log2(1.0 * transcriptomes_dict[time_t/60][gene]/transcriptomes_dict[time_0/60][gene])
return tmp
In [26]:
len(transcriptome_fold_change(time_0, time_t))
Out[26]:
In [27]:
def proteome_fold_change(time_0, time_t):
tmp = {}
for gene in results[time_0]['timecourses']:
if gene in results[time_t]['timecourses'] and results[time_0]['timecourses'][gene][-1] != 0 and \
results[time_t]['timecourses'][gene][-1] != 0:
tmp[gene] = np.log2(1.0 * results[time_t]['timecourses'][gene][-1]/results[time_0]['timecourses'][gene][-1])
return tmp
In [28]:
len(proteome_fold_change(time_0, time_t))
Out[28]:
In [29]:
def scatter_dicts_interactive(dict1, dict2, xaxislabel=None, yaxislabel=None, fontsize=None):
from bokeh.plotting import figure, show, ColumnDataSource
from bokeh.models import HoverTool
from bokeh.io import output_notebook
plotkeys = list(dict1.viewkeys() & dict2.viewkeys()) # set intersection
xs = [dict1[key] for key in plotkeys]
ys = [dict2[key] for key in plotkeys]
plt.close()
output_notebook()
source = ColumnDataSource(data=dict(x=xs, y=ys, desc=plotkeys, ))
hover = HoverTool(
tooltips=[
#("index", "$index"),
#("(x,y)", "($x, $y)"),
("", "@desc"),
]
)
p = figure(plot_width=600, plot_height=600, tools=[hover,"wheel_zoom,box_zoom,reset"],
#x_axis_type="log",
x_range=[0.8*min(xs), max(xs)/0.8],
#y_axis_type="log",
y_range=[0.8*min(ys), max(ys)/0.8],
title="",)
p.xaxis.axis_label = xaxislabel
p.yaxis.axis_label = yaxislabel
p.xaxis.axis_label_text_font_size = "10pt"
p.yaxis.axis_label_text_font_size = "10pt"
p.scatter('x', 'y', size=5, source=source, fill_color=None, line_alpha=0.5, line_color='green')
show(p)
In [30]:
times = [0, 300, 600, 900]
for i, tx in enumerate(times):
#fig, ax = plt.subplots(rows, columns, figsize=(16, 16))
for ty in times[i+1:]:
scatter_dicts_interactive(transcriptome_fold_change(tx, ty), proteome_fold_change(tx, ty),
"transcript abundance (fold change) ("+str(tx)+" s to "+str(ty)+" s)",
"proteome abundance (fold change) ("+str(tx)+" s to "+str(ty)+" s)")
5. Analysis for WCM_GRN (abundant periodic genes found by Katja)¶
In [31]:
genes_shortlist = {'HHF1': 'YBR009C', 'HHF2': 'YNL030W', 'HHT1': 'YBR010W', 'HHT2': 'YNL031C', 'HTB1': 'YDR224C',
'HTB2': 'YBL002W', 'PMA1': 'YGL008C', 'RNR1': 'YER070W', 'RNR2': 'YJL026W', 'RNR3': 'YIL066C',
'RNR4': 'YGR180C'}
In [32]:
translation_rate_shortlist = {}
for (genecommon, gene) in genes_shortlist.items():
translation_rate_shortlist[genecommon] = []
for time in np.arange(0, 3600, 300):
if gene in results[time]['timecourses']:
rate = float(results[time]['timecourses'][gene][-1]-results[time]['timecourses'][gene][-301])/300.
else:
rate = 0.
translation_rate_shortlist[genecommon].append(rate)
In [33]:
translation_ratedf = pd.DataFrame(translation_rate_shortlist, index=np.arange(0, 3600, 300))
In [34]:
pl = translation_ratedf.plot()
pl.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
pl.set_xlabel("time[s]")
pl.set_ylabel("proteins/s")
Out[34]:
In [35]:
translation_ratedf
Out[35]:
In [36]:
transcriptome_df = pd.DataFrame(transcriptomes_dict).transpose()
del transcriptome_df[0]
transcriptome_df.index = np.arange(0, 3601, 300)
transcriptome_df
Out[36]:
In [37]:
transcriptome_df_short = transcriptome_df[transcriptome_df.index<3600][genes_shortlist.values()]
In [38]:
transcriptome_df_short.columns = genes_shortlist.keys()
transcriptome_df_short
Out[38]:
In [39]:
translation_rate_per_mrna = translation_ratedf/transcriptome_df_short
In [40]:
pl = translation_rate_per_mrna.plot()
pl.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
pl.set_xlabel("time[s]")
pl.set_ylabel("proteins/mRNA/s")
Out[40]:
Translation efficiency goes down in S phase (at 1500 s) because something (ribosomes, tRNA?) becomes limiting
In [41]:
translation_rate_per_mrna
Out[41]:
In [42]:
translation_rate_per_mrna.to_csv('translation_rate_per_mrna.csv')
Comparison with initiation probabilities:translation_rate_per_mrna
In [43]:
init_rates = pkl.load(open("../parameters/init_rates_enhanced_median.p", "rb"))
init_rates_common = {gene: init_rates[genes_shortlist[gene]]*1e6 for gene in genes_shortlist}
In [44]:
scatter_dicts(dict(translation_rate_per_mrna.ix[3300]), init_rates_common, xaxislabel="translation rate per transcript",
yaxislabel="initiation rate * 1E6", lws=False, logx=False, logy=False)
Translation rate is practically proportional to initiation rate.
In [45]:
init_rates_common
Out[45]:
6. Translational efficiency over time (like 5, but at genome-scale)¶
In [46]:
translation_rate_shortlist = {}
genes = transcriptomes_dict[0].keys()[1:]
for gene in genes:
translation_rate_shortlist[gene] = []
for time in np.arange(0, 3600, 300):
if gene in results[time]['timecourses']:
rate = float(results[time]['timecourses'][gene][-1]-results[time]['timecourses'][gene][-301])/300.
else:
rate = 0.
translation_rate_shortlist[gene].append(rate)
translation_ratedf = pd.DataFrame(translation_rate_shortlist, index=np.arange(0, 3600, 300))
Translatome:
In [47]:
translation_ratedf
Out[47]:
Is the first row identically zero?
In [48]:
translation_ratedf.ix[0].describe()
Out[48]:
No it isn't. (Because the 0 refers to the interval [0 s, 300 s] etc.)
Transcriptome:
In [49]:
transcriptome_df = pd.DataFrame(transcriptomes_dict).transpose()
del transcriptome_df[0]
transcriptome_df.index = np.arange(0, 3601, 300)
transcriptome_df
Out[49]:
Divide one by the other to get the translational efficiency per mRNA:
In [50]:
translation_rate_per_mrna = translation_ratedf/transcriptome_df
translation_rate_per_mrna
Out[50]:
In [51]:
translation_rate_per_mrna.dropna(axis=1, how='all', inplace=True)
In [52]:
pl = translation_rate_per_mrna.ix[:, :18].plot()
pl.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
pl.set_xlabel("time[s]")
pl.set_ylabel("proteins/mRNA/s")
#pl.set_yscale("log")
Out[52]:
Looks like the suppressed translation efficiency is not just limited to histones.
In [53]:
translation_rate_per_mrna
Out[53]:
In [54]:
translation_rate_per_mrna.mean(axis=1, skipna=True)
Out[54]:
In [55]:
translation_rate_per_mrna.mean(axis=1, skipna=True).plot()
Out[55]:
Confirms the hypothesis that something is limiting during S phase.
7. Translational efficiency by initiation probability¶
Now we want to explode the previous graph by initiation probability.
Load unadjusted initiation probabilites:
In [56]:
init_rates = pkl.load(open("../parameters/init_rates_plotkin.p", "rb"))
Frequency distribution of initiation rates:
In [57]:
plt.hist(init_rates.values(), bins = 10**np.linspace(np.log10(min(init_rates.values())),
np.log10(max(init_rates.values())), num=50))
plt.gca().set_xscale("log")
plt.gca().set_xlim([min(init_rates.values()), max(init_rates.values())])
plt.gca().set_xlabel('initiation rates [1/s]')
plt.show()
Define translation rate as the average number of proteins produced per second over a 300 s interval:
In [58]:
translation_rate = {}
for gene in init_rates:
translation_rate[gene] = []
for time in np.arange(0, 3600, 300):
if gene in results[time]['timecourses']:
rate = float(results[time]['timecourses'][gene][-1]-results[time]['timecourses'][gene][-301])/300.
else:
rate = 0.
translation_rate[gene].append(rate)
Convert to data frame:
In [59]:
translation_rate_df = pd.DataFrame(translation_rate)
translation_rate_df.index = np.arange(0, 3600, 300)
translation_rate_df['time'] = translation_rate_df.index
translation_rate_df
Out[59]:
Convert initiation rates to data frame:
In [60]:
init_rates_df = pd.DataFrame.from_dict(init_rates, orient='index')
init_rates_df.columns = ['init_rates']
init_rates_df['init_rates'] *= 1000000
Define exponential bins for bucketing:
In [61]:
bins = 10**np.linspace(np.log10(min(init_rates_df['init_rates'])), np.log10(max(init_rates_df['init_rates'])), num=20)
bins
Out[61]:
In [62]:
translation_rate_df_long = pd.merge(init_rates_df, pd.melt(translation_rate_df, id_vars=['time']),
left_index=True, right_on='variable')
In [63]:
scatter_dicts(translation_rate_df_long['init_rates'].to_dict(), translation_rate_df_long['value'].to_dict(),
logx=False, logy=False)
Not interesting.
Cut data frame into bins:
In [64]:
translation_rate_df_long['bin'] = pd.cut(translation_rate_df_long.init_rates, bins)
translation_rate_df_long.head()
Out[64]:
In [65]:
translation_rate_df_long_grouped = translation_rate_df_long.groupby(['bin', 'time'])
translation_rate_df_long_averaged = translation_rate_df_long_grouped.mean()
Remember this is a hierarchical column index:
In [66]:
translation_rate_by_IP_by_time = translation_rate_df_long_averaged.unstack()['value']
In [67]:
xpos = np.arange(translation_rate_by_IP_by_time.shape[0])
ypos = np.arange(translation_rate_by_IP_by_time.shape[1])
yposM, xposM = np.meshgrid(ypos + 0.5, xpos + 0.5)
zpos = np.zeros(translation_rate_by_IP_by_time.shape).flatten()
In [68]:
dx = 0.5 * np.ones_like(zpos)
dy = 0.1 * np.ones_like(zpos)
dz = translation_rate_by_IP_by_time.values.ravel()
In [69]:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(12, 9))
ax = fig.add_subplot(111, projection='3d')
#values = np.linspace(0.2, 1., xposM.ravel().shape[0])
#colors = cm.rainbow(values)
ax.bar3d(xposM.ravel(), yposM.ravel(), zpos, dx, dy, dz, #color=colors,
alpha=0.5)
Out[69]:
Too ugly. We try a (first: linear) heatmap instead.
In [70]:
df = translation_rate_by_IP_by_time
plt.pcolor(df, vmin=0, vmax=2.5)
# set the limits of the plot to the limits of the data
#plt.axis([x.min(), x.max(), y.min(), y.max()])
plt.yticks(np.arange(0.5, len(df.index), 1), df.index)
plt.xticks(np.arange(0.5, len(df.columns), 1), df.columns, rotation='vertical')
plt.show()
Logarithmic heatmap:
In [71]:
df = np.log(1+translation_rate_by_IP_by_time)
plt.pcolor(df, vmin=0, vmax=2.5)
# set the limits of the plot to the limits of the data
#plt.axis([x.min(), x.max(), y.min(), y.max()])
plt.yticks(np.arange(0.5, len(df.index), 1), df.index)
plt.xticks(np.arange(0.5, len(df.columns), 1), df.columns, rotation='vertical')
plt.show()
Try again with non-exponential bins:
In [72]:
bins = np.linspace(min(init_rates_df['init_rates']), max(init_rates_df['init_rates']), num=20)
bins
Out[72]:
In [73]:
translation_rate_df_long['bin'] = pd.cut(translation_rate_df_long.init_rates, bins)
translation_rate_df_long.head()
Out[73]:
In [74]:
translation_rate_df_long_grouped = translation_rate_df_long.groupby(['bin', 'time'])
translation_rate_df_long_averaged = translation_rate_df_long_grouped.mean()
In [75]:
translation_rate_by_IP_by_time = translation_rate_df_long_averaged.unstack()['value']
In [76]:
xpos = np.arange(translation_rate_by_IP_by_time.shape[0])
ypos = np.arange(translation_rate_by_IP_by_time.shape[1])
yposM, xposM = np.meshgrid(ypos + 0.5, xpos + 0.5)
zpos = np.zeros(translation_rate_by_IP_by_time.shape).flatten()
In [77]:
df = translation_rate_by_IP_by_time
plt.pcolor(df, vmin=0, vmax=2.5)
# set the limits of the plot to the limits of the data
#plt.axis([x.min(), x.max(), y.min(), y.max()])
plt.yticks(np.arange(0.5, len(df.index), 1), df.index)
plt.xticks(np.arange(0.5, len(df.columns), 1), df.columns, rotation='vertical')
plt.show()
In [78]:
df = np.log(1+translation_rate_by_IP_by_time)
plt.pcolor(df, vmin=0, vmax=2.5)
# set the limits of the plot to the limits of the data
#plt.axis([x.min(), x.max(), y.min(), y.max()])
plt.yticks(np.arange(0.5, len(df.index), 1), df.index)
plt.xticks(np.arange(0.5, len(df.columns), 1), df.columns, rotation='vertical')
plt.show()
